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EVALUATION 

The  necessity  for  more  complex  software  systems  In  such  areas  as  command  and 
control  and  Intelligence  has  led  to  the  desire  for  better  methods  for  pre- 
dicting software  errors  and  reliability  to  insure  that  software  produced  Is 
of  higher  quality  and  of  lower  cost.  This  desire  has  been  expressed  in 
numerous  Industry  and  Government  sponsored  conferences,  as  well  as  In 
documents  such  as  the  Joint  Commanders’  Software  Reliability  Working  Group 
Report  (November  1975).  As  a result,  numerous  efforts  have  been  initiated  to 
develop  and  validate  mathematical  models  for  predicting  such  quantities  as 
the  number  of  remaining  errors  in  a software  package  and  the  time  to  achieve 
a desired  level  of  reliability.  In  addition,  efforts  have  been  initiated 
to  develop  better  methods  for  determining  when  a software  package  should  be 
released  to  a potential  user.  However,  these  efforts  have  not  produced 
measures  with  the  desired  accuracy  or  confidence  for  general  applicability. 

This  effort  was  initiated  in  response  to  this  need  for  developing  better  and 
more  accurate  software  error  prediction  and  demonstration  tests  and  fits  into 
the  goals  of  RADC  TPO  No.  5,  Software  Cost  Reduction  in  the  subthrust  of 
Software  Quality  (Software  Modeling).  This  report  summarizes  the  development 
of  mathematical  models  for  predicting  quantities  such  as  the  expected  number 
of  errors  during  both  software  development  and  software  maintenance.  These 
models  assume  errors  are  not  corrected  with  probability  1,  i.e.  imperfect 
development  and  maintenance.  The  report  also  describes  the  development  of 
statistical  tests  for  determining  whether  a software  package  should  be 
accepted  or  rejected  after  completion  of  testing.  The  importance  of  these 
developments  is  that  they  represent  the  first  attempt  to  develop  both  soft- 
ware error  prediction  models  that  incorporate  imperfect  debugging  and  thus 
more  closely  reflect  the  actual  software  error  detection  and  correction 
process,  and  software  demonstration  tests  that  allow  better  statistical 
criteria  for  accepting  a software  package. 

The  theory  and  equations  developed  under  this  effort  will  lead  to  much  needed 
predictive  measures  for  use  by  software  managers  in  more  accurately  tracking 
software  development  projects  in  terms  of  stated  error  and  reliability 
objectives.  In  addition,  the  associated  confidence  limits  and  other  related 
statistical  quantities  developed  under  this  effort  will  insure  more  wide- 
spread use  of  these  techniques.  The  acceptance  criteria  developed  will  permit 
better  control  of  the  release  of  software  packages  so  that  software  is  not 
given  to  a potential  user  before  it  is  ready  for  operational  usage.  Finally, 
the  measures  developed  under  this  effort  will  be  applicable  to  current  soft- 
ware development  projects  and  thus  help  to  produce  the  high  quality,  low  cost 
software  needed  for  today's  systems. 

ALAN  N.  SUKERT 
Project  Engineer 
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ABSTRACT 

This  report  provides  a summary  of  the  technical  activities 
pursued  under  Contract  F30602-76-C-0097  with  RADC  during  January 
1976-April  1978.  Research  work  discussed  in  previous  reports 
under  this  contract  is  summarized  and  some  additional  work  is 
described.  Also  included  is  a brief  description  of  research  in 
progress. 
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1 . INTRODUCTION 


During  the  past  ten  years  the  field  of  software  engineering 
has  grown  considerably  in  importance  and  scope.  A primary  motiva- 
tion for  this  growth  has  come  from  an  ever  increasing  cost  of 
developing  and  maintaining  software  systems.  This  is  specially 
true  for  the  DOD  which  needs  high  quality,  low  cost  software  for 
its  operations.  As  a result  of  this  increased  importance,  various 
fields  within  software  engineering  are  maturing  into  disciplines 
of  study  and  research,  for  example,  software  design  techniques, 
structured  programming  and  other  improved  programming  methodologies, 
program  testing  and  debugging  techniques,  software  performance 
modelling,  and  techniques  of  program  validation  and  verification. 

The  ultimate  objective  of  studies  in  all  these  fields  is  to  develop 
tools  that  will  be  useful  in  the  design,  development  and  operational 
phases  of  the  software  life  cycle.  The  objective  of  studies  dealing 
with  software  error  analysis  and  modelling  has  been  to  develop 
analytical  tools  which  can  be  used  for  improving  software  per- 
formance. Such  studies  can  be  classified  into  one  (or  both)  of 
two  categories.  In  the  first  category  the  emphasis  is  on  the  analysis 
of  software  error  data  collected  from  small  or  large  projects,  during 
development  and/or  operational  phases.  Studies  in  the  second 
category  are  primarily  aimed  at  the  development  of  analytical  models 
which  are  then  used  to  obtain  the  reliability  and  other  quantitative 
measures  of  software  performance. 

Typical  of  the  first  category  are  the  studies  by  Akiyama  [1) 
Belady  and  Lehman  [3],  Fries  [6],  Endres  (5],  Baker  [2],  Motley  et  al 


(18],  Miyamoto  [16],  Willman  et  al  [35],  Schneidewind  [26], 


Shooman  et  al  f 29 1 , Suker*-  [30,31),  Rye  et  al  [24],  Thayer  et  al 
[32],  and  Wagoner  (34).  These  studies  range  in  size  from  an  analysis 
of  small  data  sets  (108  errors) , e.g.  Wagoner  [34] , to  analysis  of 
large  sets  (3500  errors),  e.g.  Thayer  et  al  [32]  and  encompass  data 
from  an  on-line  system  [16],  an  operating  system  [3],  to  that  from 
the  Apollo  project  [24]  . 

In  the  second  category  of  papers,  several  models  have  been 
proposed  and  studied  during  the  last  six  years.  These  include 
'exponential  type'  models  of  Shooman  [28],  Jelinski  and  Moranda 
[11,12],  and  Schick  and  Wolverton  [25];  models  based  on  the  non- 
homogeneous  Poisson  process  proposed  by  Goel  and  Okumoto  [9]  and 
Schneidewind  [27]  , and  a Bayesian  model  by  Littlewood  and  Verrall 
[15].  Halstead  [10]  has  developed  a theory  based  on  'software  physics 
for  various  measures  of  the  performance  of  a software  system. 

Musa  [19]  has  introduced  a model  which  is  based  on  a large  number 
of  parameters  derived  from  the  software  system  being  modelled. 

Trivedi  and  Shooman  [33]  consider  a Markov  model  in  which  they 
incorporate  the  time  spent  for  removal  of  errors. 

Most  of  the  above  studies  assume  that  errors  are  removed  with 
certainty  when  detected.  The  purpose  of  the  investigation  summarized 
in  this  report  is  to  develop  and  study  models  for  software  per- 
formance which  account  for  the  probabilistic  nature  of  the  pro- 
grammer's action  during  debugging  and  operational  phases  of  the 
software  system,  to  provide  a methodology  for  classical  and  Bayesian 
inference  for  various  quantitative  measures  of  performance,  to 
develop  optimum  Bayesian  software  correctional  limit  policies,  and 
to  develop  software  reliability  demonstration  test  plans.  Results 
of  these  studies  are  summarized  in  Sections  2 through  6. 
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2.  A SOFTWARE  PERFORMANCE  MODEL  UNDER  IMPERFECT  DEBUGGING 

The  purpose  of  this  modelling  effort  was  to  establish  quanti- 
tative measures  for  software  systems  by  incorporating  the  probabilistic 
nature  of  the  programmer ' s actions  during  the  debugging  phase.  Such 
occurrences  have  been  termed  as  recurrent  errors  by  Fries  16] , 

Thayer  [32]  and  Willman  et  al  [35],  erroneous  debugging  by  Miyamoto 
[ 16 )/  bad  fixes  by  Jones  [13]  and  accounted  as  an  error  reduction 
factor  by  Musa  [19]  . In  this  study  we  call  them  the  imperfect 
debugging  errors. 

Even  though  the  presence  of  the  imperfect  debugging  phenomenon 
has  been  known,  no  published  model,  with  the  possible  exception  of 
Musa's  error  reduction  factor,  provides  an  explicit  way  to  acccunt 
for  it  in  software  performance  prediction.  The  model  and  other 
related  results  for  this  topic  are  summarized  below.  Details  of 
this  work  are  given  in  [7].  These  results  are  useful  for  software 
development  personnel  in  establishing  manpower  requirements  to 
achieve  a desired  quality  in  the  system  before  it  is  released  for 
operational  use.  Trade-off  studies  between  cost  of  debugging  and 
software  quality  can  be  undertaken  using  the  results  given  below. 

2.1  Model  and  Main  Results 

The  following  parameters  are  used  for  model  development  and 
analysis  : 

N = the  initial  number  of  errors  in  the  software  system  at 
the  beginning  of  the  debugging  activity 

p = probability  that  the  error  causing  a software  failure  is 
removed  when  detected 


4 


The  expected  time  required  to  obtain  a software  system  with 
njj  remaining  errors  is  given  by 


N-n. 


E 


>TN,n0’  ' * • 


2.1.2  Probability  distribution  of  a given  number  of  remaining 
errors  at  time  t 

Probability  that  there  are  nQ  remaining  errors  at  time  t 
is 

Pm  „ (t)  = P(X (t) =nn | X (0) =N) 


N,n 


where 


and 


= Gm  ( t)  -G  _•«  (t)  , nn-0 ,1,2, ...  ,N  , 
N,nQ  N,n0~l  0 


GN,N<tJ  1 


GN.-llt'  ' 0 ' 

Also,  the  expected  number  of  errors  reraining  at  time  t 


is 


E IX (t) | X (0) =N)  - Ne-p>t  . 


2.1.3  Expected  number  of  total  and  imperfect  debugging  errors 

The  expected  total  number  of  errors,  an<*  errors  due 

to  imperfect  debugging,  DN<t>»  during  a debugging  time  period  t 
are  given  by 


\ 

f 


V‘>  - 


3 (l-e'lpt 
P 


and 


DN(t) 


q-V^ 


2.1.4  Reliability  function 

The  software  system  reliability  between  (k-l)st  and  kth 
failures  is  given  by 

M*)  - Vik:1iPk-J-1,je-1N-|k-i-lllu  . 

K j=0  J 

2.1.5  Gamma  approximation 

The  computation  of  the  quantity  BM  . and  hence 

N • J 

G..  (t)  becomes  cumbersome  when  N is  large.  We  have  found  that 

N,n0 

the  following  gamma  approximation  yields  satisfactory  results  for 
large  software  systems  with  large  values  of  N . 


N,n, 


(t) 


‘BMBXlTi  . e-Bxax 


‘0  u r(a) 

where  the  scale  parameter  P and  the  shape  parameter  a are 
estimated  as 


6 = pX 


N 

£ 1/j 

j=n0+1 


N 

E 

j=n0+! 


1/j 


.2 


and 


N ; 

( £ 1/j) 

j=n0+1 

N 2 

£ 1/j 

j=n0+1 
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2.1.6  Numerical  examples 


a software  system  with  N=100,  A=0.02.  The  probability  distribu 


above  and  are  shown  in  Figure  2.1  for  p=0.9.  We  see  that  at  t=200 


say,  the  probability  of  having  zero  errors  in  the  system  is  approx 


imately  0.1,  of  one  error  about  0.15,  of  2 about  0.23  and  so  on 


Plots  of  expected  number  of  remaining  errors  at  various  times  are 


there  will  be  about  18  errors  left  in  the  system.  From  a study 


Figure  2.1  Probability  Distribution  of  Time 
to  i\q  Remaining  Errors 
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2.2  Analysis  of  Total  and  Imperfect  Debugging  Errors  in  a Real- 
Time  Control  System 

The  results  described  in  Section  2.1  were  used  to  analyze 
the  software  error  data  from  a large-scale  software  project  - a 
real-time  control  system  for  a land-based  radar  system  developed 
by  Raytheon  Co.  in  a modular  fashion,  see  Willman  (351.  The 
data  base  was  extracted  from  2165  Software  Problem  Reports  written 
against  109  operational  software  modules  over  the  development 
phase.  Details  of  this  analysis  are  given  in  Appendix  A. 


3.  AVAILABILITY  ANALYSIS  OF  SOFTWARE  SYSTEMS  UNDER  IMPERFECT 
MAINTENANCE 


In  this  section  we  describe  a model  developed  for  the  opera- 
tional phase  of  the  software  system  subject  to  imperfect  error 
maintenance  and  also  incorporate  the  time  spent  for  error  maintenance. 

3.1  Model  and  Performance  Measures 

The  following  parameters  are  used  in  this  model: 

N = Initial  number  of  errors  in  the  software  system  at  the 
beginning  of  software  operation 
p = Probability  that  the  error  causing  a software  failure 
is  removed/maintained  when  detected, 
q = 1-p  , the  Probability  of  imperfect  maintenance/removal 
A^  = Software  error  occurrence  rate  per  unit  time  when  there 
are  i errors  in  the  system 
= Software  error  correction/maintenance  rate  per  unit 
time  when  i errors  remain  in  the  system 

Consider  a stochastic  process  (X(t),  t>0},  whose  states  are 
defined  as 


X(t) 


i if  the  software  system  is  operational  while  there 
are  i errors  in  the  system  (i=0,l,2, . . . ,N) 

D if  the  software  system  is  down  for  error  removal 
i.e.,  for  maintenance. 

Then,  the  process  (X(t),  t>0}  forms  a semi-Markov  process  with  the 
one  step  transition  probability  (the  probability  that  the  next 
up-down  cycle,  resulting  in  j remaining  errors,  will  be  completed 
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by  time  t when  a software  package  has  i remaining  errors  at 
time  zero)  given  by 


-A^t  ->ii^ 

p (1-e  1 )•  (1-e  1 ) 

if 

j=i-l 

-X . t -p  . t 

q (1-e  1 )• (1-e  1 ) 

if 

j=i 

0 

otherwise 

Based  upon  the  above  model,  expressions  for  the  following 
performance  measures  of  the  software  system  are  derived  by  Okumoto 
and  Goel  [21] s 

(i)  Distribution  of  time  from  N to  a specified  number  of 
remaining  errors  n^ 

(ii)  Expected  time  required  for  the  system  to  go  from  N to 
nQ  errors,  EtT^l 

(iii)  Probability  of  the  system  being  operational  at  some  time 

with  n_  remaining  errors,  P„  (t)  . 

0 N,n0 

(iv)  Software  system  availability,  i.e.,  the  probability  of  the 

N 

system  being  operational  at  time  t , A(t)  = £ p„  (t) 

n0=0N'no 

(v)  Probability  that  the  number  of  errors  remaining  in  the 

system  is  n , P[N(t)=n]  and  the  expected  number  of  errors 
at  t , E [N (t) ] 

(vi)  Expected  number  of  errors  detected  and  corrected  by  t , 

D C 

denoted  by  Mjj  (t)  and  M^t)  , respectively. 

The  above  quantities  are  useful  for  software  managers  for 
estimating  the  time  and  manpower  requirements  to  achieve  a desired 
level  of  performance  during  the  operational  phase  of  the  software 
system. 
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3.2  Numerical  Example 


For  illustration  purposes  consider  the  case  when  X^=iX  , 
p^=ip,  N=100,  p=0.9,  X=0.02  and  p=0.05.  Using  the  expressions  for 
various  performance  measures  from  [21] , the  plots  of  state  occupancy 
probabilities  and  availability  are  given  in  Figure  3.1  and  the  plots 
of  and  ^(t)  are  given  in  Figure  3.2.  From  a study  of 

such  plots  for  various  values  of  X , p and  p,  one  can  obtain  adequate 
information  about  the  behavior  of  the  software  system  as  well  as 
about  the  resource  requirements  to  achieve  a desired  level  of  per- 
formance. Thus,  if  p is  known  to  be  0.9,  and  the  values  of 
X and  p that  can  be  provided  for,  then  the  system  availability 


at  t=300  will  be  approximately  0.5.  If  this  is  not  satisfactory, 


one  has  to  provide  additional  or  better  resources  that  will  yield 
better  values  for  one  or  more  parameters. 


A Nomogram  for  the  Expected  Time  to  a Specified  Number  of 
Errors  and  to  Determine  Manpower  Requirements 


We  present  a simple  nomogram  to  calculate  the  expected  time 
required  to  remove  a specified  number  of  errors  from  a software 
system  which  will  satisfy  the  desired  performance  requirements. 

We  consider  the  case  when  X.=iX  and  p.=ip,  i.e.,  when  the  error 


of  remaining  errors  with  X and  p , respectively,  the  constants 
of  proportionality.  Letting  the  ratio  X/p*p,  we  have 


NO.  OF  DETECTED  8 CORRECTED  ERRORS 


N =100 
X = 0.02 

fi-  0.05 
p = 0.9 


A nomogram  which  gives  the  contours  of  4>(<(>,p)  in  the  (*,p) 


To  illustrate  the  use  of  this  nomo 


consider  the  case  when  N=100 


per  day  and  the  desired  value  of  n 


100+1 

10+1 


Step  2.  Corresponding  to  <(>=.963  and  p=0.4,  from  Figure  3.3  the 


value  of  *(.963,0.4)  is  3.2 


Step  3.  Compute  the  value  E ( T 


Thus,  for  the  given  conditions  the  expected  time  required 
to  go  from  100  to  10  errors  is  177.8  days. 


Now  we  consider  another  application  of  this  nomogram  to 


show  how  it  can  be  used  in  determining  the  manpower  requirements 
for  a specified  objective  of  remaining  errors.  Suppose  we  want 
to  go  from  N=1000  to  nQ=10  errors  in  E(T1000  ^Q)=500  days  when  the 
error  occurrence  rate  is  X=0.01  errors  per  day  and  the  probability  p 


the  manpower  requirement  to  accomplish  this  objective. 

The  first  thing  to  determine  is  the  value  of  p that  will 


We  know  that 


* (<KP)  = pXE[TM  ] 


and  hence 


*U,p)  = (.95)  (.01)  (750)  = 7.125. 

From  the  nomogram  in  Figure  3.3,  for  * ( 4> , p ) =7 . 125  and 
(t»=log(^Yji)=1.96,  the  value  of  p=0.58  and  hence  p=X/p=0 . 0172 
errors/day.  If  the  error  removal  rate  per  person  per  day  is  0.01, 
then  we  will  need  1.7  people  to  meet  the  desired  objective. 
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4.  BAYESIAN  AND  CLASSICAL  INFERENCE  FOR  THE  IMPERFECT  DEBUGGING 
AND  MAINTENANCE  MODELS 


In  this  section  we  describe  two  methods  for  statistical  infer 
ence  of  the  parameters  of  the  models  described  in  Sections  2 and  3. 
The  first  one  is  the  classical  approach  based  on  maximum  likelihood 
estimation  and  the  second  is  a Bayesian  approach  based  on  the  prior 
distributions  of  the  unknown  parameters.  The  parameters  under 
consideration  are  the  initial  number  of  software  errors  N,  the 
error  occurrence  rate  for  each  error  X , and  the  probability  of 
perfect  debugging  p . An  additional  parameter  for  the  imperfect 
maintenance  model  is  y . We  give  only  the  method  for  the  model 
of  Section  2.  The  same  procedure  can  be  used  for  the  model  of 
Section  3. 

The  available  data  for  estimation  purposes  is  generally 
given  as  t= (t^, t2 , . . . , t^) , the  times  between  software  failures 
and  £-(y^,y2> • • • *y n>  > an  indicator  variable  for  imperfect  debugging 
such  that  y^=l  if  the  ith  failure  is  caused  by  an  error  due  to 
imperfect  debugging  and  y^=0  if  the  error  is  not  due  to  imperfec 
debugging. 

The  maximum  likelihood  and  Bayesian  approaches  to  the  estima- 
tion of  the  above  parameters  are  summarized  below.  The  details 
of  the  procedure  are  given  in  [20] . 


4.1  Maximum  Likelihood  Method 

The  likelihood  function  for  N,  p and  X is 

-{N-p(i-l) )Xt 


L(N,p,X I t,y)=  n {N-p (i-1) }e 
i=l 


yi  1-yi 

i , q(i-l)  , 1 . N- (i-1) i 

‘{N-^<r-lT}  *{N-F(i-ir 


} 

l 

- 

* 


1 
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The  maximum  likelihood  estimates  (N,p  and  X)  are  obtained  as  solu- 
tions of  the  simultaneous  non-linear  equations 

n n 1-y . 

X l t.  = E jT-TT1^ 
i-1  1 i-1  N‘(l“1} 

n n 

X E (i-l)t.  = E y./q 
i=l  1 i=l  1 

n 

n/X  = E {N-p (i-1) } t . . 
i=l  1 

The  joint  100(l-a)%  confidence  regions  for  N,  p and  X are 
obtained  from 

1 2 


i (N,p,X | t,£)-i (N,p,X | t,£)  = ^ x3;a  » 


where 


*(• |t,y)  = logL ( • |t;£)  . 

The  estimated  variance-covariance  matrix  for  N,  p and  X is 

-1 


rNN 

rNp 

rNX  " 

E 

cov 

rpN 

rpP 

rPx 

- rXN 

rxP 

rxx  - 

N=ft 

p=p 

X = X 


where 


rNN  * t 1/{N- (i-1) ) {N-p (i-1) ) 


IT  a r a A 

rNp  pN 
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rNA  = rAN  = A l/tN-p(i-l)> 

1 n 

i £ (i-1) /{N-p (i-1) } if  q^O 
q i=l 

« if  q=0  . 

rpx  - rxp  - - r .^li-D/tK-Pd-l)) 

rAX  ■ "''l2  • 

4.2  Bayesian  Inference 

Now  we  describe  a Bayesian  approach  for  obtaining  posterioj 
point  estimates  and  the  highest  posterior  density  (HPD)  region  fc 
parameters  N,  p and  A . 

The  choice  of  the  prior  distribution  for  a parameter  is 
governed  by  several  factors.  In  our  case  we  take  the  conjugate 
priors,  which  for  N and  A are  gamma  distributions  while  for  j 
it  is  a beta  distribution,  i.e.. 


P(N)*  N°-1  e“eN  , 

N>0 

P(p)-  pff-1(l-p)p_1  , 

0<p<l 

P(A)«*  A11-1  e“YX  , 

A>0  . 

By  applying  Bayes  theorem  the  joint  posterior  distributic 
of  N,  p and  A for  given  priors  and  the  data  is  obtained  as 


P (N,p,  A | t,j£)o».  p(N,p,A)L<N,p,A|t,£) 


Let  N,  p and  A be  the  Bayesian  point  estimates  for  N,  p 
and  A , respectively.  That  is,  the  point  (N,p,A)  is  the  mode  of 
the  joint  posterior  distribution  p(N,p, A | t,yj , i.e.  it  attains 
its  maximum  at  (N,p,A).  Then,  N,  p and  A are  the  values  that 


satisfy 


n n 1_yi  a-l 

"X  iii  ^ + ill  + “S"  ‘ 6 = ° ' 


A r^(i-l)ti-Eyi/(l-p)  + ^ = 0 * 


n/A  - z {N-p (i-1) } t . + - y = 0 . 

i=l  1 A 

Finally,  the  100(l-a)%  Bayesian  confidence  region  is  given 


f(N,p,A)  = C , 


where 


f(N,p,A)  = nlogA  - Z {N-p(i-l)}t. 

i=l  1 


n n 

+ Z y.log(l-p)  + r (1-y. ) log{N- (i-1) } 
i=l  1 i=l  1 


+ (a-l)logN-eN 


+ (y-DlogA-yA 


+ (n-l)logp  + (p-1) log (1-p) 


C » f(N,p,A)  - | Xj. B 
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5.  BAYESIAN  SOFTWARE  CORRECTION  LIMIT  POLICIES 

The  objective  of  the  investigation  was  to  provide  an  optimur 
correction  limit  policy  for  a large-scale  software  system  subject 
to  random  error  occurrences  and  error  removals  in  an  operational 
phase.  When  an  error  occurs  a corrective  action  is  undertaken  to 
remove  it.  Such  an  action  can  be  scheduled  at  two  levels,  which 
we  call  Phase  I and  Phase  II.  By  Phase  I we  mean  that  the  cor- 
rective action  will  be  undertaken  by  the  programmer  while  Phase  II 
action  is  undertaken  by  a system  analyst  or  system  designer.  First, 
Phase  I corrective  action  is  scheduled  for  a specified  time  T . 

If  the  error  is  not  corrected  in  this  time,  it  is  referred  to 
Phase  II.  This  sequence  of  corrective  actions  in  an  operational 
phase  is  shown  in  Figure  5.1.  Our  objective  is  to  determine  the 
optimum  value  T*  of  T which  minimizes  the  long  run  average 
cost.  Two  models  are  developed  for  this  purpose.  In  the  first 
model  we  assume  that  the  cost  of  observations  of  error  occurrence 
and  correction  time,  prior  to  the  implementation  of  the  optimum 
policy,  is  negligible.  The  second  model  incorporates  the  cost  of 
observations. 


6.  SOFTWARE  RELIABILITY  DEMONSTRATION  TEST  PLANS 

The  purpose  of  this  task  was  to  describe  the  theory, 
methodology,  and  procedures  for  software  reliability  demonstration 
tests.  Such  tests  are  to  be  conducted  to  ensure  that  the  software 
system  being  considered  for  acquisition  has  achieved  desired  reli- 
ability. The  situation  we  have  in  mind  for  applying  these  tests 
is  that  the  software  system  has  gone  through  the  usual  development 
phases,  including  testing  and  debugging,  and  is  being  presented 
for  testing  to  demonstrate  its  reliability.  Software  reliability 
for  this  purpose  will  be  defined  as  the  probability  that  a given 
software  program  operates  for  a given  time  period,  without  a soft- 
ware error,  on  the  machine  for  which  it  was  designed  given  that 
it  is  used  within  design  limits. 

A complete  description  of  the  results  on  this  project  is 


given  in  Appendix  B. 
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APPENDIX  A 


AN  ANALYSIS  OF  RECURRENT  SOFTWARE  ERRORS 
IN  A REAL-TIME  CONTROL  SYSTEM 

In  this  Appendix  we  present  an  analysis  of  software  error 
data  from  a large-scale  software  project  using  the  imperfect 
debugging  model  discussed  in  Section  2.  The  model  parameters 
are  estimated  from  the  data  and  the  values  predicted  from  the  model 
are  compared  with  the  observed  values.  Joint  confidence  regions 
for  the  parameters  are  also  constructed  which  permit  a study  of 
the  sensitivity  of  predictions. 

A. 1 Analysis  of  Error  Data  from  a Real  Time  Control  System 

A real-time  control  system  for  a land-based  radar  system 
was  developed  by  Raytheon  Co.  in  a modular  fashion.  Nearly  all 
of  the  modules  were  written  in  J0VIAL/J3 . The  error  data  base 
was  extracted  from  2165  Software  Problem  Reports  (SPRs)  written 
against  109  operational  software  modules  over  the  development 
phases  and  is  described  in  [35].  Table  A.l  shows  the  distribution 
of  the  SPRs  by  month  opened  during  a 22  month  period  of  integration, 
acceptance  and  operational  testing  phases. 

The  available  data  give  the  number  of  total  errors  (u^)  and 
the  number  of  imperfect  debugging  or  recurrent  errors  (v^)  detected 
by  time  t^  . The  parameters  under  consideration  are  the  initial 
number  of  software  errors  (N) , the  error  occurrence  rate  for  each 
error  (X) , and  the  probability  (p)  of  perfect  debugging.  We  first 
estimate  the  parameters  N,  p and  X from  these  data  (t,u ,v)  and 
then  compare  the  results  obtained  from  I DM  with  the  observed  values. 
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Using  the  data  in  Table  A.l  and  che  results  from  Section  2 


the  results  are  obtained  as  shown  in  Table  A. 2 


The  joint  confidence  regions  for  N and  X for  p=0.974  are 


by  month  are  shown  in  Figure  A. 2,  and  the  plots  of  actual  and 


predicted  number  of  remaining  errors  are  given  in  Figure  A. 3 
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APPENDIX  B 

SOFTWARE  RELIABILITY  DEMONSTRATION  TEST  PLANS 

BL1.  INTRODUCTION  AND  PROBLEM  DEFINITION 

The  purpose  of  this  report  is  to  describe  the  theory,  methodology, 
and  procedures  for  software  reliability  demonstration  tests.  Such 
tests  are  to  be  conducted  to  ensure  that  the  software  system  being 
considered  for  acquisition  has  achieved  desired  reliability.  The 
situation  we  have  in  mind  for  applying  these  tests  is  that  the 
software  system  has  gone  through  the  usual  development  phases,  in- 
cluding testing  and  debugging,  and  is  being  presented  for  testing 
to  demonstrate  its  reliability.  Software  reliability  for  this 
purpose  will  be  defined  as  the  probability  that  a given  software 
program  operates  for  a given  time  period,  without  a software  error, 
on  the  machine  for  which  it  was  designed,  given  that  it  is  used 
within  design  limits. 

In  order  to  develop  an  appropriate  test  plan,  we  must  first 
choose  a model  which  adequately  describes  the  error  occurrence 
phenomenon.  Most  software  error  models  are  based  on  the  assump- 
tion that  successive  errors  follow  a decreasing  failure  rate  be- 
cause of  the  reduction  in  the  number  of  remaining  errors.  However,  in 
this  appendix  we  assume  that  the  times  between  errors  during  the 
demonstration  phase  follow  an  exponential  distribution,  l.e.  have 
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a constant  failure  rate.  As  a justification  for  this  assumption, 
we  contend  that  the  demonstration  time  and  the  number  of  errors  en- 
countered during  demonstration  will  be  relatively  small  and  hence  a 
constant  failure  rate  model  will  be  an  adequate  representation  of 
the  error  phenomenon.  At  worst,  this  model  will  yield  a somewhat 
conservative  test  from  the  viewpoint  of  DOD. 

Let  the  distribution  of  error  occurrence  times  be  given  by 

f (t | A>  - A exp  (-tA),  t >_  0,  A > 0 (B.l-1) 

Then  the  following  measures  of  software  reliability  can  be  used 
interchangeably : 

, Reliability,  R(t)  =>  P(t  > t)  = e_xX  (B.l-2) 

. Meantime  Between  Software  Failure  (MTBSF)  * 1/A 

. Software  Failure  Rate  (SFR)  = A 

In  the  following  we  shall  take  A as  the  measure  of  software 
reliability. 

In  the  classical  set  up,  the  demonstration  tests  are  generally 
designed  to  distinguish  between  two  values  of  A,  namely  the  maximum 
acceptable  SFR,  A^,  and  the  specified  SFR,  A^.  The  decision  to 
accept  or  reject  the  software  system  is  based  upon  test  results 
which  are  subject  to  random  fluctuation.  A loss  is  incurred  when 
either  an  accept  or  a reject  decision  is  wrongly  taken. 

The  loss  corresponding  to  wrong  decisions  is  quantified  by  two 
risks  called  the  producer's  risked  - p(Rl\)),and  the  consumer's  risk 
6 ■ P(A|A1),  where  A and  R denote  acceptance  and  rejection  of  the 
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software  system,  respectively 


In  this  context  various  types  of  demonstration  plans  can  be 
designed  to  limit  a and  6 to  desired  values.  For  example,  a 
truncated  single  sample  plan  for  the  system  is  employed  as  follows. 

The  plan  consists  of  using  the  software  in  an  environment  which  is 
representative  of  the  operational  environment  for  a time  period  T. 

The  number  of  errors  r encountered  during  this  time  period  is 
recorded.  (In  this  study  we  assume  that  all  errors  are  of  the  same 
severity.  When  errors  are  classified  according  to  the  degree  of 
severity,  the  problem  becomes  quite  complicated  and  is  beyond  the 
scope  of  this  investigation).  If  r is  less  than  or  equal  to  a 
pre-specified  number,  r*,  the  software  is  accepted.  Otherwise,  the 
system  is  rejected.  The  design  of  such  a plan  consists  of  obtaining 
the  quantities  T and  r*  such  that  the  desired  risks  a and  3 are  satis- 
fied. 

One  disadvantage  in  the  above  approach  is  that  X is  assumed 
to  be  an  unknown,  fixed  constant,  In  practice  there  may  be  suffi- 
cient reason  to  believe  that  knowledge  about  \ is  available  in  some 
quantifiable  form  and  one  would  like  to  Incorporate  such  informa- 
tion in  the  development  of  the  demonstration  test  plan.  Another  impor- 
tant situation  arises  when  the  risks  associated  with  the  demonstration 
test  are  not  adequately  represented  by  the  classical  risks  (atP)  end 
Interest  lies  in  associating  risks  with  the  posterior  distribution  of 
In  such  cases  one  resorts  to  a Bayesian  approach  for  test  design. 
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The  problem  under  consderatlon  Chen,  is  the  design  of  single 
sample  software  reliability  demonstration  plans  when  error  occur- 


I 


rence  times  are  assumed  to  follow  the  exponential  distribution. 

Before  describing  the  development  of  the  test  plans,  we  first 

l 

discuss  the  various  risks  that  arise  in  the  above  situations  in 
Section  B.2.  Fixed  time,  classical  tests  are  then  considered  in 
Section  B.3.  In  Sections  B.4  and  B.5  Bayesian  tests  are  developed 
for  two  situations:  (i)  X has  a noninforma tive  prior  distribution, 
and  (ii)  information  about  X can  be  quantified  from  the  testing 
and  debugging  data. 


Section  B.6  presents  the  step-by-step  procedure  to  be  employed 
to  demonstrate  the  properties  and  performance  of  the  demonstration 
test  plans. 


B..  2.  DEFINITIONS  AND  INTERPRETATION  OF  RISKS 


The  following  quantities  are  of  interest 

1.  The  probability  P(r|X  ” Xq)  - a that  a software  systea  with 
specified  SFR  is  rejected. 

2.  The  probability  P(a|X  - X^)  * 6 that  a software  system  with 
maximum  acceptable  SFR  is  accepted. 

3.  The  probability  P(a|X  X^)  * B that  a software  system  which  is  of 
unacceptable  reliability  is  accepted. 

4.  The  probability  P (R | X X^)  » a that  a software  system  of  acceptable 


reliability  is  rejected. 

5.  The  probability  P(X  ^ X^|a)  ■ 3*  that  the  SFR  of  a system  which  has 
been  accepted  is  more  than  the  maximum  acceptable  SFR. 

6.  The  probability  P(X  _>  Xq|a)  - 3**  that  the  SFR  of  a system  which 
has  been  accepted  is  more  than  the  specified  SFR. 

7.  The  probability  P(X  _<  Xq|r)  « a*  that  the  SFR  of  a system  which 
has  been  rejected  is  less  than  the  specified  SFR. 

8.  The  probability  P(R)  that  the  software  system  is  rejected. 

9.  The  probability  P(A)  that  the  software  system  is  accepted. 

10.  The  probability  P(X  >_  X')  that  a-priori,  the  software  system 

has  SFR  which  is  core  than  X'. 

These  quantities  will  now  be  discussed  in  some  detail. 
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B.  2.1,  Classical  Risks  (a, 3) 


The  classical  producer's  risk  a and  consumer's  risk  B are  defined 
as  follows: 


{ 

i 


o “ P(r|A  - Aq),  the  probability  of  rejecting  a (B.2-1) 

software  system  whose  SFR  is  equal  to* 
the  specified  value,  A . 

B = P(A|A  * A^) , the  probability  of  accepting  a software  (B.2-2) 
system  whose  SFR  Is  equal  to  the  maximum 
acceptable  value. 

The  (a, 3)  risks  represent  two  points  on  the  classical  operating 
characteristic  (OC)  curve  which  Is  a plot  of  P(A|A)  versus  A.  These 
risks  do  not  provide  an  explicit  control  of  the  probability  of  accep- 
tance for  values  of  A other  than  A^  and  Aq.  However,  P(a|A)  decreases 
monotonically  with -A1.  Hence,  if  A < Aq,  the  probability  of  rejection 
is  less  than  a,  If  A > A^,  the  probability  of  acceptance  is  less 
than  3.  The  shape  of  the  OC  curve  governs  the  degree  of  protection 
provided  in  the  indifference  zone  between  A^  and  Aq. 

B.  2.2.  Average  Risks  (a,B) 


The  average  risks  are  defined  as  follows 

a » P(R|A  Aq),  the  probability  of  rejecting  a software  (B.2-3) 
system  with  a SFR  less  than  or  equal  to 
the  specified  value, Aq. 

B ■ P(A|A  A^),  the  probability  of  accepting  a software  (B.2-4) 
system  with  a SFR  greater  than  or  equal 
to  A^. 
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Mathematically,  the  risks  may  be  expressed  as: 


P(R|X  < Aq) 


P(R,A  < Aq) 

P (A  < A0) 


fU 

I P(Rj A)p(A)dA 


O 


p(A)dA 


(B. 2-5) 


a =»  J P(R|A)*p(A|A  < AQ)dA  (B.2-6) 

0 

and,  similarly 

00 

0-J  P(A|A).p(A|A  > Ax)dA  (B.2-7) 

x1 

The  average  risks  provide  the  following  protection.  If  the 
producer  produces  a large  number  of  software  systems,  then,  in  the 
long  run,  less  than  a percent  of  the  desired  ones  will  be  rejected. 

If  the  consumer  buys  a large  number  of  software  systems,  then,  in 
the  long  run,  less  than  100  8 percent  of  the  bad  systems  will  be 
accepted.  No  explicit  control  on  the  probability  of  acceptance  is 
provided  at  any  specific  value  of  A when  we  use  the  average  risk 
criteria. 


I 


2 . 3 Posterior  Risks  (a*,  g*) 

The  (a*,  P*)  risks  are  defined  as  follows 

a*  - P(A  £ Aq|r) 


This  risk  is  the  long  run  probability  of  a rejected  software  system 


being  good. 


p*  = p(a  > ax |A) 


This  risk  is  the  long  run  probability  of  an  accepted  system  being  bad. 


Mathematically , 


a*  = P(A  < Aq|R)  = | P'(A|R)dA  - 


V 

| P(R|A) 


p(A)dA 


|°° P(R|  A)p(A)dA 


OD 

f P(A|A)  P(A)dA 

CO  J 

f X1 

e*  - P(A  > AX|A)  - P(A|A)dA  - — 

^1  P(A|A)  p(A)dA 

■ 0 

These  risks  can  also  be  Interpreted  in  a "degree  of  belief" 
sense.  Thus , a*  would  represent  a persons's  degree  of  belief  that 
if  a software  system  has  been  rejected,  it  has  a SFR  which  is  better 
than  the  specified  value.  Similarly,  6*  would  be  the  degree  of  belief 
that  the  SFR  of  an  accepted  system  is  worse  than  the  maximum  acceptable 
value. 


B-8 


B»  2.4  Probability  of  Rejection  P(R) 

This  is  a single  number  given  by 

OO 

P(R)  - j P(R|A)p(A)dA 
0 


(B.2-12) 


or 


P (R)  - 1 - P (A) 


,oo 


1 - 


P(A| A)p(A)dA 


(B.2-13) 


Note  that  the  integration  is  over  the  entire  range  of  A and  specifi- 
cation of  A0> which  is  usually  specified  in  conjunction  with  a producer’s 
risk,  is  unnecessary. 

In  the  frequency  sense  we  have 


„ Total  number  of  systems  rejected 
Total  number  of  systems  tested 


For  the  producer  this  criterion  implies  that,  in  the  long  run,  less 
than  (100) *P(R)  percent  of  the  systems  will  be  rejected. 


B.2.5.  Alternate  Posterior  Consumer's  Risk  R** 
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B.  3.  DESIGN  OF  TEST  PLAN  FOR  CLASSICAL  RISKS 


Now  we  consider  the  design  of  a single  sample  plan  where  the 
parameters  are  T and  r*.  Since  the  time  to  software  failure  is 
exponential,  the  observed  number  of  software  errors  r in  fixed  time 
T has  a Poisson  distribution,  i.e. 


f (r  | A)  =*  - 


-TA  (TA)' 


r! 


(B.3-1) 


The  two  risks  can  be  written  as: 

-TA 


r*  e (TA  )r 

y o_ 

n r! 

r=0 


- 1 - a 


(B.3-2) 


and 


•TX1 

r*  e (TA,)r 


I T!-5- 

r-0  r! 


- 6 


(B»  3-3) 


Given  A^,  Aq,  a and  $,  the  above  equations  can  be  simultaneously 
solved  to  obtain  software  test  time  T and  allowable  number  of  errors 
r*. 

Individual  specification  of  A^  and  AQ  is  not  necessary.  Let 
K “ A1/Aq,  and  let  T*  ■ T Aq.  Then,  to  satisfy  the  stated  risks  we 
can  write 


I 
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M - a 


(8,3-4) 


r! 


and 


r*  -KT* 

l- 


r-0 


(KT*) 

r! 


r 


l e 


(B.3-5) 


Given  (K,  a,  $)  the  above  two  equations  can  be  solved  to  obtain 
(T*,  r*) . These  equations  can  be  solved  numerically  or  by  using 
tables  of  cumulative  Poisson  probabilities.  Clearly,  test  plans 
with  Identical  A^A^  have  the  same  T*  and  r*  values.  Given  the 
specified  value  AQ,  the  actual  test  time  is  obtained  as  T - t*Aq* 
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4.  BAYESIAN  TEST  PLANS  FOR  NONINFORMATIVE  PRIOR 
DISTRIBUTION 

In  this  section  we  consider  the  case  when  it  is  possible  to 
express  the  unknown  parameter  A in  terms  of  a metric  <f>(A)  so  that 
the  corresponding  likelihood  is  data  translated.  This  means  that 
the  likelihood  function  for  $(A)  is  completely  determined  a-priori 
except  for  its  location  which  depends  on  the  software  failure  data 
yet  to  be  observed.  This  state  of  indifference  can  be  expressed 
by  taking  <P(\)  to  be  locally  uniform,  and  the  resulting  prior  dis- 
tribution is  called  noninformat ive  for  <p(l)  with  resepect  to  data. 

A more  detailed  discussion  on  noninformat ive  prior  distribution 
tan  be  found  in  Box  and  Tiao  [1]*  . In  our  case,  a noninformat ive 

prior  for  X is 

p(X)  « A~1/2  (B.4-1) 

Now,  if  T is  the  test  time  for  a software  system,  the  number  of 
errors  in  T will  be  given  by  a Poisson  distribution  with  parameter 
TX  as  mentioned  earlier.  Letting  A - TX,  the  prior  for  A can  be 
written  as 

P(A)  « A"1/2  (B.4-2) 


* References  at  the  end  of  this  Appendix 


The  joint  probability  distribution  of  A and  r is: 


L 


p(r,A)  - p(r |A) *p(A) 


(_A)r  • e~A  , -1/2 

r!  A 


or 


p(r,A)  = b 


, A 


r-l/2  . -A 


r! 


(B,4-3) 


where  b'  is  the  normalizing  constant. 

To  get  the  marginal  distribution  of  r,  we  let  A be  defined  over  the 
range  (0,d)  where  d is  sufficiently  large  and  d < °°.  Then  we  have 


0° 

p(r)  * j P(r,A/*dA  = j°V- 
0 


, r-l/2  . -A 


r! 


dA 


(BA-4) 


Now  we  choose  d such  that 


lo 


, .r-l/2  . -A 


b’  A 


e_dA  <e 


r ! 


(B A-5) 


for  some  given  sufficiently  small  G > 0:’  Then 


P(r)  - T (r  +\) 


(B.4r6) 


and 


P(A) 


r* 

l P(r) 

r-0 


b*r(r  + 4 ) 


l 

r-0 


r! 


(Bi4-7) 
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Iii  order  to  get  expressions  for  various  risks,  we  proceed 


as  follows. 

ft 

Let  Aq  = TXq.  Then  from  (B.4-3)  we  get, 

? 

^q*  accept)  “ P(A  A^,  accept) 

t 

i 

Substituting  the  above  expressions  into  the  appropriate 
formulae  in  Section  2,  we  get  the  equations  for  desired  risks. 

For  example,  if  we  are  interested  in  the  risk  combination  (a,fJ*),  we 
get  from  (B.2-3)  and  (B.2-9),  respectively: 

5*1-  P(A|X  < 10) 

P(A,  A<\0) 

= 1 

P(A  < A ) 


3*  * P(X  > X1|A) 


1 - p(X  < X |A) 


r*  r1 

l ) b*Ar 

r-0  L 


-1/2  -A. 

• e dA 


r*  b'  r (r  + f) 
r-0  71 


r*  ;A 

I f A" 

r-n  J 


-1/2  -A. 

• e dA 


3*  = 1 - 


r*  T(r  + 

r — - 

r-0  r! 


(a4-10) 


If  the  consumer's  risk  of  Interest  Is  3**,  then  from  (B.2-14) 


we  have 


i 


3**  - P(X  > Xq|a) 

- 1 - P(X  < X_|A) 


r* 

l { Ar 

r-0  { 


~1/2  * e‘A  dA 


(B.  4-11) 


r*  r (r+fr 

r-0  r! 

B-  16 


The  design  plan  T and  r*  is  obtained  by  simultaneously  solving 
two  equations,  one  each  for  the  producer's  and  the  consumer’s  risk. 
The  proper  choice  of  risk  combinations  will  be  governed  by  the 
protection  desired  and  agreed  upon  by  the  vendor  and  the  buyer  of 
the  software  system  or  systems. 

Two  nomograms  for  the  design  of  plans  for  (5,$*)  and  are 

given  on  the  following  pages, 
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FIG.  B-l  COUNTOURS  OF«.£*FOR  DESIGN  OF  TEST  PLANS 


FIG.  B-2  CONTOURS  OF a,/9"*FOR  DESIGN  OF  TEST  PLANS 


B.  5.  BAYESIAN  TEST  PLANS  FOR  CONJUGATE  PRIOR 


In  some  situations,  data  on  failures  and  times  between  failures 
duringJthe  testing  and  debugging  phase  may  be  available.  Such  data 
can  be  used  to  obtain  an  appropriate  prior  distribution  for  X,  the 
software  failure  rate.  A flexible  prior  and  one  which  is  mathe- 
matically tractable  in  this  case  is  a two  parameter  gamma  given  by 


TP  .pT  -tA 
p(X)  =*  pp  * x • e 


(B.5-1) 


Estimates  of  the  parameters  p and  t can  be  obtained  from  the  available 
data.  Using  this  prior,  the  expressions  for  the  various  producer's 
and  consumer's  risks  are  obtained  as  follows: 


Expressions  for  Producer's  Risks: 


Classical: 


Average : 


r*  e TA°  • (TA  )r 

1 - “ “ l 7j  ~ 

r-0  r* 

"I  e‘T*  • (TA)2 


(B.5-2) 


1 - a 


L e • 

r-0  . t"  ,%p-l  . -tA 

' ~ r.  W X C " 

X0 

| ^ •X'’'1  ' ''TX'dX 


(B.5-3) 
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Posterior: 


Probability  of  Rejection: 


P(R) 


Expressions  for  Consumer's  Risks : 


Classical: 


e-  l 

r=0 


-TKA 

r*  e u • (KTA0)r 


r! 


Average : 


(B.  5-5) 


(B.  5-6) 


dA 


(B.  5-7) 


I 
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B.  6.  PROCEDURE  FOR  DEMONSTRATING  THE  PERFORMANCE 
OF  THE  TEST  PLANS 

This  section  provides  a step-by-step  procedure  that  will  be  used 
to  demonstrate  the  properties  and  performance  of  the  demonstration  test 
plans  for  a software  system  as  developed  in  the  previous  sections. 

These  procedures  will  be  followed  in  conducting  the  demonstration  on 
site. 

Prior  to  conducting  the  demonstration  test  in  practice,  the  DOD 
must  make  sure  that  ' he  software  system  has  been  developed  according 
to  specifications  and  all  the  stated  development  phases  have  been 
carried  out  to  meet  the  objectives. 

The  steps  to  be  followed  for  a demonstration  test  are  as  follows. 

1.  Decide  the  risk  criteria  (producer's  and  consumer's  risks)  for 
demonstration. 

2.  Choose  the  values  of  the  two  risks. 

3.  Design  a test  plan  (T,r*)  which  meets  the  risks  in  Step  2 as 
closely  as  possible.  (Designs  will  be  obtained  using  computer 
programs  developed  for  this  purpose.) 

4.  Simulate  software  errors  using  the  software  error  simulator. 

5.  Apply  the  demonstration  test  of  Step  (3)  to  simulated  errors 
and  make  accept/reject  decisions. 

6.  Compare  the  recorded  results  with  theoretical  results. 


B-  23 


B.?.  REFERENCES 


Box,  G.E.P.  and  Tiao,  G.C.  (1973) 

Bayesian  Inference  in  Statistical  Analysis. 
Goel,  A.L.  and  Joglekar,  A.M.  (1976) 


Addison-Wesley 


"Reliability  Acceptance  Sampling  Plans 
Distribution:  Risk  Criteria  and  Their 
RADC-TR-76-294,  Volume  II.  (A033516) 


Based  Upon  Prior 
Interpretation, " 


APPENDIX  C 


COMPUTER  PROGRAMS 

C.l  Programs  for  the  Imperfect  Debugging  Model  (Section  2) 

Computer  programs  to  compute  the  following  quantities,  for 

given  N,  p and  A are  given  in  Tables  C.l  to  C.8. 

• Mean  and  variance  of  the  first  passage  time  from  N to  n^ 

(SUBROUTINE  FRSTRT) 

• Pdf  and  cdf  of  the  first  passage  time  from  N to  n^ 
(SUBROUTINE  FPT2) 

• Probability  of  the  software  system  having  n^  errors  remaining 
at  time  t (SUBROUTINE  STATE) 

• Expected  number  of  errors  remaining  at  time  t and  expected 
number  of  total  errors  and  imperfect  debugging  errors  detected 
by  time  t (SUBROUTINE  MEAN) 

• MTTF  and  reliability  for  given  k at  time  x 
(SUBROUTINE  TBF) 

The  programs  are  self-explanatory  and  include  the  required 

subroutines . 
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SIS  P19S  ZS  BEST  QUALITY  PRACTICABLE 

IM  OOPY  fUUUSHSD  TO  DOG  - 


TABLE  C.l  SUBROUTINE  FRSTPT 


subroutine  frstpt( 
function  - 


usage 

parameters 


p - 

r - 
nO  - 
et  - 
vt  - 


n»p»  r »nO»et»vt»a»b) 

compute  mean  and  variance  of  first  passage  time 
from  n to  nO  t and  estimate  shape  and  scale 
parameters  of  Gamma  distribution 
call  frstpt  (n*p»rrnOfet»vtFa»b) 

(input.)  initial  no.  of  errors 

(input.)  prob.  of  perfect  debugging 

(input.)  detection  rate  / error 

(input.)  desired  no.  of  errors 

(output.)  mean  time  from  n to  nO 

(output.)  variance  of  time  from  n to  nO 

(output.)  shape  parameter  of  Gamma  distribution 

(output.)  scale  parameter  of  Gamma  distribution 


subroutine  frstpt 
xl«0.0 
x2=0.0 
nl=nO+l 
do  55  J=nl »n 
xl=xl+l . 0/f loat( J) 
x2=x2+1.0/float(J> 
55  continue 
et=xl/p/r 
vt=x2/p/p/r/r 
b=et/vt 
a=et*b 
return 
>end 


(n»Pfr»nO»et»vt»a»b) 


i 


TABLE  C.2 


SUBROUTINE  FPT2 


c 

subroutine 

fpt2 

c 

c 

function 

c 

usage 

c 

parameters 

n 

c 

p 

c 

r 

c 

nO 

c 

t 

c 

pdf 

c 

cdf 

of  first  passage  tine  from 


- compute  p.d.f.  and  c.d.f< 
n to  nO 

- call  fpt2  <n»p» r ?nO»t»pdf »cdf ) 

- (input.)  initial  no.  of  errors 

- (input.)  prob.  of  perfect  debugging 

- (input.)  detection  rate  / error 

- (input.)  desired  no.  of  errors 

- (input.)  tine 

- (output.)  p.d.f.  of  first  passage  tine  from 

- (output.)  c.d.f.  of  first  passage  tine  fron 
c read,  subroutines  - frstpt»  mdgann 


to  nO 
to  nO 


c 

c 


subroutine  fpt2  (niPi r»nO»t»pdf *cdf ) 
call  frstpt  (mprr»nO»et»vt»a»b) 
x=b*t 

xl=(a-l .0)*alog(x) 
x2=al gamma (a) 
x3=xl-x-x2 

if  (x3  .It.  -88.0)  go  to  33 
pdf=exp(x3> 
go  to  44 
33  pdf =0.0 

44  call  ndganm  (x.arcdf) 
return 
22  cdf=0.0 
pdf=0.0 
return 
end 
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TABLE  C. 3 


SUBROUTINE  STATE 


subroutine  state  (n.p. r» nO.t.pt)- 


e 

c 

c 

c 

c 

e 

c 

c 

e 

c 


function 

usage 

parameters 


- compute  probability  of  being  nO  errors  remaining 

at  time  t 

- call  state  (n.p. r »nO. t.pt) 

n - (input. ) initial  no.  of  errors 
p - (input.)  prob.  of  perfect  debugging 
r - (input.)  detection  rate  / error 
nO  - (input.)  desired  no.  of  errors 
t - (input.)  time 

pt  - (output. ) prob.  of  being  nO  errors  at  time  t 
• fpt2 


read,  subroutines 


subroutine  state  (n»p» r»nO* t.pt) 
if  (nO.eo.n)  go  to  33 
call  fpt2(n»p»r»n0»t»pdf »cdf ) 
ptl=cdf 
go  to  44 
33  ptl=l .0 
44  nl=n0-l 

if  (nl)  11.22.22 
11  Pt 2-0.0 
go  to  99 

22  call  fpt2  (n.p. r.nl .t.pdf »cdf ) 
pt2=cdf 
99  pt=ptl-pt2 
return 
end 


b 
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TABLE  C . 5 SUBROUTINE  TBF 


subroutine  tbf  (nipi r r k » x r rel » xmttf ) — 

fund i on  - compute  reliability  and  mttf 

usade  - coll  tbf  (n»p* r»k»xi rel .xmttf ) 

parameters  n - (input* > initial  no.  of  errors 

p - (input.)  prob,  of  perfect  debuddind 
r - (input.)  dec  tion  rate  / error 
k - (input.)  k-th  failure 
x - (input. > time 

rel  — (output.)  reliability  at  time  x after  (k-l)st  failure 
xmttf  - (output.)  mean  time  between  (k-l)st  and 

k-th  failures 


subroutine  tbf  (n»p» r»k»x» rel .xmttf ) 

xmttf=l .0/(f loat(n)-p*f loat(k-l) )/r 

rel=exp(-x/xmttf ) 

return 

end 
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TABLE  C. 6 SUBROUTINE  MDGAMM 


subroutine  mddamm  (x.p.prob)" 


-from  imsl- 


e function  - compute  incomplete  damma  distribution 

C usade  - call  nddamm  (x>PrProb) 

C parameters  x -(input.)  value  to  which  damma  is  to  be  intedrated 

C p — (input.)  damma  parameter 

prob  - (output.)  prob.=intedral  of  damma(p)  to  x 
read,  subroutines  - damma 


subroutine  mddamm  (x.p.prob) 
dimension  v(&)ivl(6) 
eouivalence  (v(3)»vl(l)) 
prob-0.0 

if  (x  .de.  0.0)  do  to  5 
do  to  9000 

5 if  (p  .dt.  0.0)  do  to  10 
do  to  9000 

10  if  (x  .ea.  0.0)  do  to  9005 
fnld=aldamma(p) 
cnt=p*alod(x> 
xcnt=»x+fnld 

if  ((cnt-wcnt)  .dt.  -88.0)  do  to  15 
ax*0 . 0 
do  to  20 

15  ax=exp(cnt-xcnt) 

20  bid>l.e35 
cut=l .e-8 

if  ((x  «le.  1.0)  .or.  <x  .It.  p))  do  to  40 

v*l .0-p 

z-x+x+1 .0 

cnt=0.0 

v(l>=1.0 

v(2)=x 

v(3)=x+l .0 

v(4)=z*x 

prob=v(3)/v<4) 

25  cnt=cnt+1.0 
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wcnt=y*cnt 

v(5)=vl ( 1 )*z-v( 1 ) Cvcnt 

v(A)=vl (2)*z-v(2)*vcnt 

If  (v(6)  .ea.  0.0)  do  to  50 

ratio=v(5)/v(6) 

reduc=abs( prob- ratio) 

if  (reduc  .dt.  cut)  do  to  30 

if  (reduc  >le.  ratiotcut)  do  to  35 
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TABLE  C.6 


(Continued) 


30  prob*ratio 
ao  to  SO 

35  prob*l .0-probtax 
bo  to  9005 
40  ration 
cnt=>l  .0 
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prob»1.0 

45  ratio*ratio+l .0 
cnt=cnt*x/ ratio 
prob=prob+cnt 

If  (cnt  .at.  cut)  ao  to  45 
prob=prob*ax/p 


50 


ao  to  9005 
do  55  i=ii4 
v(i>*vl(i) 


55  continue 

if  (abs(v(5)>  .It.  bia)  ao  to  2 
do  A0  i=l »4 


v(i>sv(i)/bia 
60  continue 
ao  to  25 


9000  continue 
9005  return 
end 


function  alaamma(p) 
if  <e  .at.  31.0)  ao  to  15 
call  aamma  (nap) 
alasniina=aloa<  a*) 
return 

15  zl=(e-0.5)*aloa<p)-p+0.5*aloa<2.0*3.1415) 
z2*l. 0/12. 0/p 
z3=l • 0/360 • 0/p/p/p 
z4=l .0/1260 • 0/p/p/p/p/p 
z5=l • 0/1680 . 0/p/p/p/p/p/p/p 
alda»ea>zl-fz2-z3+z4-z5 
return 
end 
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TABLE  C.7  SUBROUTINE  GAMMA 


subroutine  M»u  (xx*ax)  fro*  imsl- 

function  - coarute  a danaa  function  of  parameter  xx 

usaae  - call  Sanaa  (xxiSx) 

parameters  xx  - (input* ) parameter  of  Sanaa  function 

Sx  - (output*)  value  of  Sanaa  function 


subroutine  sanna(xxiSx) 
if(xx-57.)  6*6*4 
4 Sx*1.0e30 
return 
6 x»xx 

err*>l  *e-4 
sx*l. 

if (x-2. ) 50*50.15 
10  if  (x-2.)  110*110*15 
15  x**x— 1 • 

SX«SX*X 
So  to  10 

50  if  (x-l.)  60*120*110 

60  if  (x-err)  62*62.80 
62  x*float(int(x))-x 

if  (abs(v)-err)  120*120*64 
64  if  (1. -x-err)  120*120*70 

70  if (x-l.)  80*80*110 
80  Sx>Sx/x 
x*x+l. 

So  to  70 
110  W“x-1 . 

av-1 . ♦ y* ( -0 .5771017+x* ( 0 . 9858540+y* ( -0 . 8764218+x* ( 0 . 8328212+y* ( -0 . 5684729 
\c+x* (0.2548205+y* (-0.051 4993 ) > ) ) > ) ) 

SXBSXtSU 
120  return 


i 


C.2  Programs  for  Simulation  of  Imperfect  Debugging  Data 

Computer  programs  to  simulate  the  data  required  to  estimate 

the  parameters  N,  p and  A , are  given  in  Tables  C.8  to  C.ll. 
These  programs  perform  the  following  functions: 

• Simulate  data  (t,j£)  for  given  N,  p and  A (SUBROUTINE  SMLT) 

. Compute  the  mle ' s of  N,  p and  A given  the  data  (t,£)  and 

also  obtain  the  estimate  of  variance-covariance  matrix 
(SUBROUTINE  MLE) 

• Compute  the  Bayesian  estimates  of  N,  p and  A for  given 
data  (t,^)  (SUBROUTINE  BAYES) 

The  programs  are  self-explanatory. 


TABLE  C. 8 


SUBROUTINE  SMLT 


c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  smlt 
function  - 

usage  - 

parameters  n - 

p - 
r - 
nn  - 
iseed  - 


t - 
lv  - 


(n.p. r.nn. iseed. t» iy) 

simulate  time  between  s/w  failures  for 
imperfect  debugging  model 
call  smlt  (niP»ri nn » i seed »tf is) 
(input.)  initial  no.  of  s/w  errors 
<input.)  prob.  of  perfect  debugging 


(input.)  detection  rate  / error 

(input.)  no.  of  observations  for  s/w  failure  time 
(input.)  an  integer  value  in  the  exclusive 
range  ( 1 . 2147483647) . iseed  is  replaced  by 
a new  iseed  to  be  used  in  subseauent  calls, 
output  vector  of  length  nn.  containing  time 
between  s/w  failures 

output  vector  of  length  nn.  indicating  the 


t'ype'  of 'error' ‘which  is  1 if  the'i-th 
failure  is  caused  by  an  error  due  to 
imperfect  debugging » or  0 otherwith. 


c 

c 

c 


read,  subroutine  - SSub 


subroutine  smlt  (nrPt r .nn. iseed.t. iy) 

dimension  t (nn) » iy (nn) . r r (2) 

nr=n 


ie-0 

do  5 i=l.nn 

if  (nr  .ea.  0)  go  to  99 
xm=l .0/r/f 1 oat (nr ) 
call  ggub  ( iseed. l.rr) 
t(i  )«=-xm#alos(  rr(  1 ) ) 
call  ggub  (iseed. l.rr) 
if  (rr(l)-(l.O-p))  22.22.33 
33  nr=nr-l 
go  to  44 
22  ie=ie+l 

44  call  ggub  (iseed. l.rr) 

if  (rr(l)-float(ie)/float(r.r)) 
66  iy( i )=0 
go  to  77 
55  iy( i >=1 
ie=ie-l 

77  print  100. i »t( i ) . iy ( i ) .nr . ie 
100  format( iS.f 15.5.315) 

5 continue 
99  return 
end 
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TABLE  C . 9 SUBROUTINE  MLE 


I 

f 

I 


t- 


e 

c 

c 

e 

c 

e 

c 

c 

c 

c 

c 

c 

c 

e 

e 

c 

c 

c 

c 


subroutine  ale  <t?iy?nn»en»ep?er»eeov) 

function  - estimate  unknown  parameters  n>  pi  and  lambda 

and  also  estimate  variance-covariance 
matrix  for  idm 

usage  — call  mle  (t?  iamnieniepi er? ecov) 

parameters  t - ( input •>  a vector  of  length  nn?  containing  time 

between  s/w  failures 

iv  - (input.)  a vector  of  length  nnr  indicating  the 
twpe  of  error  which  is  1 if  the  i-th 
failure  is  caused  by  an  error  due  to 
imperfect  debugging?  or  0 otherwith 


n — (input.)  no.  of  observations  for  s/w  failure  time 
en  — (output.)  an  estimate  of  parameter  n 

ep  - (output.)  an  estimate  of  parameter  p 

er  - (output.)  an  estimate  of  parameter  lambda 

ecov  - (output.)  an  estimate  of  variance-covariance 

matrix  (3x3) 


subroutine  mle  <t?iu?nn?en?ep?er?ecov) 
dimension  t <nn) ? ia(nn) ?ecov(3?3) 
xl-0.0 
x2=0.0 
v=0.0 

do  5 i=l?nn 
xl=xl+t< i ) 

x2=x2+t<i)*float<i-l> 
float ( iu < i ) ) 

5 continue 

enO=float(nn+l> 
epO= 1 .0 
do  IS  J=1 ?20 
x3=xl*enO-epO*x2 
x4«=0.0 
x5=0.0 
do  25  i=l?nn 
x4=x4+f loat < 1-iy < i ) )/<en0-f loat( i-1 ) )**2 
x5=x5+f loat( l-iy( i ) )/(en0-f loat(i-l ) ) 

25  continue 

f0=x5-xl*float(nn)/x3 

fn®-x4+xl*xl*f loat(nn)/x3/x3 

f p*-xl*x2#f loat <nn)/x3/x3 

ph0=f loat(nn)*< 1 .0-ep0) *x2/x3-g 

rhn*=-<  1 .0-ep0)  *f loat  (nn)  4x24x1  /x3/x3 

Php«float(nn)*x2*(-1.0+(1.0-ep0)*x2/x3)/x3 

hk“fn*php-phn*fp 

hh»- ( f 0*php-ph0*f p )/hk 

xk»- ( f n*ph0-phn*f 0 ) /hk 

en0=en0+hh 

ep©=epO+xk 


.Ctv 


s* 
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TABLE  C.9 


(Continued) 


print  100» Jf enO.epO 
100  for»at<i5'?el5.5> 

if  (aeaxl <sbs<hh/enO) t abs<xk/epO) ) .le.  0.00001)  do  to  11 
15  continue 
11  en*enO 
er=eeO 

•r*f loat(nn)/x3 
rnn=0 . 0 
rpr=0.0 
rnr=0.0 
do  3S  i = l »nn 

rnn=rnn+l .0/(en-f loat< i-1 ) )/(en-ep*f loat( i-1 ) ) 
rpr=rpr+f loat< i-1 >/<en-ep#f loatt i-1 ) ) 
rnr=rnr+l .0/<en-ep*f loat( i-1 > ) 

35  continue 
rnr=rnr/er 
rpp=rpr/< 1 .0-ep) 
rpr=-rpr/er 
rrr=f loat(nn)/er/er 

rx=rnn*rpp*rrr-rnr*rpp*rnr-rnn*M>r*rpr 
ecov(l 1 1 ) = < rpptrri — rprtrpr >/rx 


ecov( 1 * 2)=rpr*rnr/rx 

ecov< 1 »3)=-rpp*rnr/rx 

ecov(2»2)=(rnn*rpr-rnr*rnr)/rx 


ecov(2»3)=-rnn#rpr/rx 

ecov(3»3)=rnn*rpp/rx 

ecov(2> 1 )=ecov< 1 p 2) 

ecov(3r 1 )=ecov( 1 »3> 

ecov(3*2)=ecov(2p3) 

return 

end 


I 
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TABLE  C.10  SUBROUTINE  BAYES 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  bases  ( t»  is»nn» alpha fbeta» pi  » rho  »>:muF3amma»  bn.  bp.br) 

function  - obtain  basesian  estimates  of  unknown  parameters  n»  p» 

and  lambda  for  idm 

usa3e  - call  bases  < t . is.nn.alpha .beta .pi  » rho.xmu. aamma .bn.bp  »br  ) 

parameters  t - < input.)  a vector  of  lenath  nn.  containina  time 
between  s/w  failures 

is  - (input.)  a vector  of  lenath  nn.  indicatina  the 
tspe  of  error  which  is  1 if  the  i-th 


nn 

alpha 

beta 

Pi 

rho 

xmu 

aamma 

bn 

bp 

br 


failure  is  caused  bs  an  error  due  to 
imperfect  debuaaina.  or  0 otherwith 
(input.)  no.  of  observations  for  s/w  failure  time 
(input.)  shape  parameter  of  aamma  prior  for  n 

(input.)  scale  parameter  of  aamma  prior  for  n 

(input.)  first  parameter  of  beta  prior  for  p 
(input.)  second  parameter  of  beta  prior  for  p 
(input.)  shape  parameter  of  aamma  prior  for  lambda 

(input.)  scale  parameter  of  aamma  prior  for  lambda 

(output.)  basesian  estimate  of  n 

(output.)  basesian  estimate  of  p 

(output.)  basesian  estimate  of  lambda 


subroutine  bases  ( t » is» nn* alpha* beta iPi » rho . xmu. aammarbnrbp.br ) 

dimension  t(nn).is(nn) 

xl«=0.0 

x2=0.0 

s=0.0 

do  5 i=l»nn 
xl*=xl+t(  i ) 

x2=x2+t(i >*f loat( i-1 ) 
wswtf loat ( is< i ) > 

5 continue 

bnO*float(nn+l ) 

bp 0*1 .0 

do  IS  J»1.20 

x3*xl*bn0-bp0*x2+damma 

X4-0.0 

xS-0.0 


014 


TABLE  C.10 


(Continued) 


do  25  i=l»nn 

x4=x4+f loat< l-i«<i ) >/<bnO-f lost ( i-1 > )**2 
x5=x5+f loat< l-iu< i ) )/<bnO-f loat< i-1 > > 

25  continue 

f0=x5-xl*< xmu+f loat<nn-l ) )/x3+<alpha-l .0>/bn0-beta 
fn=-x4+xl*xl*<  xmu+f 1 oat < nn-1 > >/x3/x3-(alrha-l  .0)/bn0/bn0 
fp=-x4*x2*(xniu+f  loat(nn-l  ) )/x3/x3 

phO=  < xmu+f loat ( nn- 1 ) >*< 1 .0-bp0)*x2/x3-y+<pi-l .0>*< 1 .0-bp0)/bp0-<  rho-1 .0) 
phn=-(l  ,0-bp0)*<xmu+f loat<nn-l ) )*>:2#xl/x3/x3 

php=<xmu+f loat(nn-l ) >*x2*<-l .0+< 1 ,0-bp0)*x2/x3)/x3-<pi-l .0)/bp0/bp0 
hk=fn*php-phn*f p 
hh=-  ( f 0*php-ph0*f  p ) /hk 
xk=-<fn*phO-phn*f 0)/hk 
bnO=bnO+hh 
bpO=bpO+xk 
print  100 » J»bnO » bpO 
100  format( i5*2el5.5) 

if  (ama>:l(abs(hh/bn0) »abs(xk/bp0> ) .le.  0.00001)  do  to  11 
15  continue 
11  bn=bn0 
bp=bp0 

br=< xmu+f loat(nn-l > )/x3 

return 

end 


jfilS  PAGE  IS  BEST  QUALITY  PRACTICABL! 
PftOS  OQPY  TURBUSHBD  TO  DDQ 


ms  Wfli  IS  B*ST  QUALITY  mCttCAM*! 
jy  acm  y\BBttSHBD  10  " 


TABLE  C.ll  SUBROUTINE  GGUB 


c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 


subroutine  ggub 
function  - 

usage  - 

parameters  iseed  - 


n - 
r<n>  - 


(iseedmtr) from  imsl 

basic  uniform  <0rl>  pseudo-random  number 
Generator 

call  ggub  (iseedmrr) 

(input.)  an  integer  value  in  the  exclusive 
range  < 1 »2147483647> . iseed  is  replaced  by 
a new  iseed  to  be  used  in  subseouent  calls, 
(input.)  no.  of  deviates  to  be  generated 
(output  vector  of  length  n»  containing  the 
deviates  in  (0»1> 


subroutine  ggub(  iseedmi  r ) 
dimension  r(l) 
double  precision  zrdpmrdp 
drm=dble< float ( 2**31 -1 ) ) 
dr*dble<f loat< 1/2**31 ) ) 
z- iseed 
do  5 i=l»n 

Z»dmod< 16807 . d0*z . dem ) 

5 r(i)=z/dpm 
iseed~z 
return 
end 


flMMNN 
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C.3  Programs  for  the  Imperfect  Maintenance  Model  of  Section  3 

Computer  programs  to  compute  the  following  quantities  of 
interest  for  given  N,  p,  X and  y are  given  in  Tables  C.12  to 
C.18. 

• Mean  and  variance  of  the  first  passage  time  from  N to  nQ 
(SUBROUTINE  COMP) 

• Pdf  and  cdf  of  the  first  passage  time  from  N to  nQ 
(SUBROUTINE  FIRST) 

• Probability  of  the  software  system  being  operational  with  nQ 
remaining  errors  at  time  t (SUBROUTINE  STT) 

• Software  system  availability  at  time  t (SUBROUTINE  AVAIL) 

• Expected  number  of  errors  detected  and  corrected  by  time  t 
(SUBROUTINE  EXPCT) 

The  programs  are  self-explanatory  and  include  the' required 
subroutines  (MDGAMM  and  GAMMA) . 


TABLE  C.12 


SUBROUTINE  COMP 


c 

subroutine 

comp  (r 

c 

function 

c 

c 

c 

usade 

c 

parameters 

n 

c 

p 

c 

r 

c 

xm 

c 

kl 

c 

k2 

c 

et 

c 

vt 

c 

a 

c 

b 

e 

rl 

c 

r2 

compute  mean  and  variance  of  first  passage 
tine?  estimate  the  gamma  parameters*  and 
obtain  the  constants  rl  and  r2 
call  comp  <n* p» r »xm*kl »k2*et * vt» a»br rl » r2) 
(input.)  initial  no.  of  errors 
(input.)  prob.  of  perfect  debuddind 
(input.)  detection  rate  / error 
(input.)  correction  rate  / error 
(input.)  first  destination 
(input.)  second  destination 
(output.)  mean  first  passage  time 
(output.)  variance  of  first  passage  time 
(output.)  shape  parameter  of  gamma  distribution 
(output.)  scale  parameter  of  gamma  distribution 
(output.)  smaller  constant 
(output.)  larder  constant 


subroutine  comp  (n»p* r *xm»kl *k2»et»vt*a»br rl*r2) 

kkl=kl+l 

kk2=k2+l 

rr=sort<  <r+xm)**2-4.0*p*r*xm) 

rl=( r+xm-rr)/2.0 

r2=( r+xm+rr )/2.0 

erl=0.0 

er2=0.0 

vrl=0.0 

vr2=0.0 

if  (kkl  .dt.  n)  do  to  11 
do  5 i=kkl»n 
erl=erl+l .0/f loat< i ) 
vrl*vrl+l .0/f loat< i >/f loat ( i ) 

5 continue 

11  erl*=erl/rl  , 

vrl=vrl/rl/rl  rtj/k 

if  (kk2  .dt.  n)  so  to  22 
do  15  i*kk2*n 
er2=er2+l. 0/f  loat(i) 
vr2"vr2+l .0/f loat ( i ) /f loat ( i > 

15  continue 
22  er2«er2/r2 
vr2*>vr2/r2/r2 
et«erl+er2 
vt“vrlkvr2 

if  (vt  . ea.  0.0)  do  to  33 
b*et/vt 
So  to  99 
33  b=0.0 
99  a=ct#b 
return 


TABLE  C. 13 


SUBROUTINE  FIRST 


subroutine  first  (n»p»r»xm»kl»k2»t*pdf »cdf »rl»r2> 


c subrout 
c function 
c 

c usage 
c parameters 


parameters  n 

p 

r 

xm 

kl 
k2 
pdf 
P cdf 
rl 
r2 

reod.  subroutines 


- compute  pdf  and  cdf  of  first  passage  time 

and  also  obtain  the  constants  rl  and  r2 

- call  first  <n?p» r »xmf kl »k2»t »pdf tcdf r rl t r2> 

- (input.)  initial  no.  of  errors 

- (input.)  prob.  of  perfect  debugging 

- (input.)  detection  rate  / error 

- (input.)  correction  rate  / error 

- (input.)  first  destination 

- (input.)  second  destination 

- (output.)  Pdf  of  first  passage  time 

- (output.)  cdf  of  first  passage  time 

- (output.)  smaller  constant 

- (output.)  larger  constant 

- comp . mdgammr  gamma 


subroutine  first  (n» p» r t xm»kl rk2» t . pdf rcdf t rl . r2> 
if  (minO(kl »k2)  .It.  0)  go  to  22 
call  comp  (n»p» r »xm»kl »k2»et * vt » ar b» r 1 » r2) 
if  (b  .eo.  0.0)  ao  to  11 
x=b*t 

xl=(a-l . 0>*alos(>:> 

K2=alsamma(a) 
x3=xl-x-x2 

if  (k3  .It.  -88. 0)  go  to  33 
pdf=e>.'P(x3) 

go  to  44 
33  pdf =0.0 

44  retirnd3amm  <x,3,cdf>  1HIS  PACB  IS  BEST  QUALITY PRACIICABL! 

11  cdf=l  .0  TR0I  OOPY  TURSISHBD  TO  DOC  ^ — 

pdf=l .0 
return 
22  cdf =0.0 
pdf =0 . 0 
return 
end 
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TABLE  C.14 


SUBROUTINE  STT 


subroutine  stt  (n»r»r»x»»nO»t»rrob>- 


e 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 


function 

usage 

parameters 


n 

p 

r 

x* 

nO 

t 

prob 


read,  subroutine 


- compute  the  probability  of  being  nO  s/u 

errors  remaining  at  time  t 

- call  stt  <n»p» r »xm» nOrtrProb) 

- (input.)  initial  no.  of  errors 

- (input.)  prob.  of  perfect  debugging 

- (input.)  detection  rate  / error 

- (input.)  correction  rate  / error 

- (input.)  specified  no.  of  errors 

- (input.)  time 

- (output.)  prob.  of  being  nO  errors  remaining 

at  time  t 

- first 


TABLE  C. 15 


SUBROUTINE  AVAIL 


subroutine  avail  (n»r» r at) 

function  - compute  s/w  system  availability  at  time  t 

usage  - call  avail  (n*p»  r»xm»  t»at) 

parameters  n - (input*)  initial  no.  of  errors 

p - (input.)  prob.  of  perfect  debugging 
r - (input.)  detection  rate  / error 
xm  - (input. ) correction  rate  / error 
t - (input.)  time 

at  - (output.)  availability  at  time  t 
read,  subroutine  - stt 


subroutine  avail  (nrpirixmit.at) 

at *0.0 

nl«n+l 

do  5 i=l*nl 

n0=i-l 

call  stt  (n»p»r»>:m»nOrt»prob) 
at=at+prob 
5 continue 
return 
end 


IRXI  FAGHt  IS  BEST  QUALITY  FRACTICABL| 
lt*&*  QOPY  finilSHS)  TO  DDQ 


IBIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 

nm  oan  puwish sd  to  dog 


TABLE  C.16  SUBROUTINE  EXPCT 


e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  expct 
function 

usase 

parameters  n 

p 

r 

x» 

t 

xmd 

xac 


read,  subroutine 


(nrP»rr>:mrt>  xmd » xmc  > 

- compute  mean  no.  of  errors  detected  and 

corrected  by  time  t 

- call  expct  (n»p» r »xm»t»xmd»xmc) 

- (input.)  initial  no.  of  errors 

- (input.)  prob.  of  perfect  debugging 

- (input.)  detection  rate  / error 

- (input.)  correction  rate  / error 

- (input.)  time 

- (output.)  mean  no.  of  errors  detected  by 

time  t 

- (output.)  mean  no.  of  errors  corrected  by 

time  t 

- first 


subroutine 

hl=0.0 

h2=0.0 

h=0.0 

do  5 i = l*n 
k=i-l 

call  first 
hl=hl+cdf 
call  first 
h2=n2+cdf 
call  first 
h=h+cdf 


expct  (n»p»r»xm»t»xmd»xmc> 


(n»p»r»xm»k»i ft » pdf »cdf » rl » r2) 
(nrPf  r»xniri»K»ti pdf  »cdf  » r 1 » r2 ) 
(n » p ? r » xm  r k > k * t f pdf  rcdf  t rl  r r2 ) 


S continue 
xmc=h/p 

xmd= ( hi  * ( 1 . 0-xm/ r 1 ) +h2*  < xm/ r2- 1 . 0 ) ) *r / ( r 1 - r2  > 

return 

end 
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TABLE  C . 17  SUBROUTINE  MDGAMM 


c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  mdsamm  (xinp rob ) from  imsl 

function  - compute  imcomplete  aamma  distribution 

usaSe  - call  mdsamm  (xiPfprob) 

parameters  x - (input*)  value  to  which  samma  is  to  be  integrated 

p - (input*)  samma  parameter 
prob  - (output.)  prob.  =intearal  of  3em»a(p)  to  x 
read,  subroutines  - samma 


subroutine  mdsamm  (x»p»prob) 
dimension  v(6)rvl(6) 
eouivalence  (v(3)rvl(l)) 
prob=0.0 

if  (x  .se.  0.0)  so  to  5 
So  to  9000 

5 if  (p  .St.  0.0)  So  to  10 
So  to  9000 

10  if  (x  .eo.  0.0)  so  to  9005 
fnls=alsamma(p) 
cnt=p*alos(x) 
wcnt=x+fnls 

if  ((cnt-ycnt)  .St.  -88.0)  So  to  15 

ax=0.0 

So  to  20 

15  ax=exp( cnt-ycnt) 

20  bis=l . e35 
cut= 1 . e-8 

if  <<x  .le.  1.0)  .or.  (x  .It.  p>)  So  to  40 

y=l . 0-p 

z=x+y+l .0 

cnt=0.0 

v(l)=1.0 

v(2)=x 

v(3)=x+l .0 

v(4)«z*x 

prob=v(3)/v<  4) 

25  cnt=cnt+1.0 
v«=y+l  .0 
z-z+2.0 
ycnt=v*cnt 

v(5)=vl(l )*z-v< 1 )*ycnt 
v(A)=vl (2)*z-v(2)*ycnt 
If  (v(6)  .eo.  0.0)  So  to  SO 
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TABLE  C. 17 


(Continued) 


ratiosv(5)/v(6) 
reduc=abs(prob-ratio) 
if  (reduc  .at.  cut)  ao  to  30 
if  (reduc  . le.  ratiotcut)  ao  to  35 
30  prob=ratio 
ao  to  50 

35  prob=l .0-prob*ax 
ao  to  7005 
40  ratio=p 
cnt=l .0 
prob=l .0 

45  ratio=ratio+l  .0 
cnt=cnt*x/ratio 
prob=prob+cnt 

if  (cnt  .at.  cut)  ao  to  45 
prob=p  rob*ax/p 
ao  to  9005 
50  do  55  i=l »4 
v(i)=vl<i) 

55  continue 

if  (abs(v(5) ) .It.  bia)  ao  to  25 
do  60  i = l r 4 
v(i)=v(i )/bia 
60  continue 
ao  to  25 
9000  continue 
9005  return 
end 


| 


function  alaamnia(p) 
if  (p  .at.  31.0)  ao  to  15 
call  aarnma  (p.ap) 
elaaMaBaloa(ap) 
return 

15  zl=<p-0.5)*aloa<e)-p+0.5«aloa<2.0*3.1415) 
z2=l .0/12.0/p 
z3=l .0/360. 0/p/p/p 
z4=l. 0/1260. 0/p/p/p/p/p 
z5*l . 0/1680 . 0/p/p/p/p/p/p/p 
alda*tika=zl+z2-z3+z4-z5 
return 
end 


fttis  FAQB  ZS  BIST  QUALITY  PRACTICABU 
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TABLE  C. 18  SUBROUTINE  GAMMA 


subroutine  damma  (xxijx) fro*  imsl 

function  - compute  a damma  function  of  parameter  xx 

usade  - call  damma  (xxrdx) 

parameters  xx  - (input.)  parameter  of  damma  function 

dx  - (output.)  value  of  damma  function 


subroutine  damma(xx.dx) 
if (xx-57. ) 4.4.4 
4 dx=l . 0e30 
return 
6 x=xx 

err=l .e-4 
dx»l  • 

if (x-2. ) 50.50.15 
10  if  (x-2.)  110.110.15 
15  x«*x-l. 

dx»dx*x 
do  to  10 

50  if  (x-1.)  40.120.110 

60  if  (x-err)  42.42.80 
42  w=float( int(x) >-x 

if  (abs(v)-err)  120.120.44 
44  if  (l.-v-err)  120.120.70 

70  if (x-1.)  80.80.110 
80  dx*dx/x 
x«x+l . 
do  to  70 
110  v*x-l. 

dy-1 . +y* ( -0 . 5771017+y* ( 0 . 9858540+y* ( -0 . 874421 8+y* (0 . 8328212+y* ( -0.5484729 
\cfw* ( 0 . 2548205+y* ( -0 . 0514993 ))))))> 
dx*dx*dv 
120  return 
end 
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C . 4 Program  for  Bayesian  Software  Correction  Limit  Policies 
(Section  5) 

Computer  programs  to  compute  the  optimum  policies  for  model  1 
and  model  2 are  given  in  Tables  C.19  to  C.23.  These  programs  per- 
form the  following  functions: 

• For  model  1,  simulate  error  occurrence  time  and  correction  time 
and  then  compute  Bayesian  estimates  of  mean  correction  time, 
optimum  correction  limit  time,  and  its  minimum  cost  per  unit 
time  (SUBROUTINE  MDLl) 

• For  model  2,  simulate  error  occurrence  time  and  correction  time 
and  then  compute  Bayesian  estimate  of  mean  correction  time, 
optimum  ocrrection  limit  time,  and  its  minimum  cost  per  unit 
times,  and  also  provide  the  optimum  sample  size 

(SUBROUTINE  MDL2) 

The  programs  are  self-explanatory  and  include  the  required 
subroutines  (eg.  DATA  1,  GGUB  and  OPTMM) . 
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ASSESS®®0 


TABLE  C. 19 


SUBROUTINE  MDLl 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

e 

c 

c 

c 

c 

c 

c 

c 


subroutine  mdll 
function  - 


usage  - 

parameters  iseed  - 


nn  - 
r - 
xmul  - 
xmu2  - 
c - 


alpha  - 


beta  - 


xn  - 
*n  - 

tt  - 
ct  - 

read*  subroutines 


(iseed*  nn  * r * xmu  If  xmu2  r c * a 1 pha  * be  t a * xn  * un  * 1 1 * c t ) 

simulate  error  occurrence  time  and  correction  time 
and  then  compute  bauesian  estimate  of  mean 
correction  time*  optimum  correction  limit  time* 
and  its  minimum  cost  per  unit  time 
call  mdll  ( iseedtnn*  r* xmul *xmu2*c* alpha  * beta r xn* an* tt*ct) 
(input.)  an  integer  value  in  the  exclusive  range 
(1 >2147483647) . iseed  is  replaced  by  a new  iseed 
to  be  used  in  subsecruent  calls. 

(input.)  no.  of  observations 
(input.)  error  occurrence  rate 
(input.)  mean  correction  time  in  phase  1 
(input.)  mean  correction  time  in  phase  2 
(input.)  a vector  of  length  3 

c(l>  - cost  per  unit  time  of  error  correction 
in  pase  1 

c(2>  - cost  per  unit  time  of  error  correction 
in  pase  2 

c(3>  - sampling  cost  per  sample  size 
(input.)  a vector  of  length  2 

alpha(l)  - shape  parameter  of  inverted  gamma 


prior  for  mean  correction  time 
in  pase  1 

alpha(2)  - shape  parameter  of  inverted  gamma 
prior  for  error  occurrence  rate 
(input.)  a vector  of  length  2 

beta(l)  - scale  parameter  of  inverted  gamma 
prior  for  mean  correction  time 
in  pase  1 

beta(2)  - scale  parameter  of  inverted  gamma 
prior  for  error  occurrence  rate 
(output.)  bayesian  estimate  of  error  occurence  rate 
(output.)  bauesian  estimate  of  mean  correction 
time  in  pase  1 

(output.)  optimum  correction  limit  time 
(output.)  minimum  cost  per  unit  time 
- datalf  optmm 


subroutine  mdll  ( iseed* nn>  r > xmul *xmu2*c*  alpha *beta*xn*yn*tt*ct> 

dimension  c(3>  * alpha (2) *beta<2) *xy<2) 

call  datal  ( iseed*nn* r *xmul *xx *yy ) 

xn*=(xx+beta<2)  )/( alpha(2)+f  loat  (nn-1 ) ) 

a <*0.0 

b=xn 

call  optmm  (nn*yy>c*alpha>beta»xmu2*a*b*yn*tt*ct) 

return 

end 
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TABLE  C. 20  SUBROUTINE  MDL2 


e 

e 

c 

c 

e 

c 

e 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

e 

c 

c 

c 

c 

c 

e 

e 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  ndl2  ( iseed. r > xnul>xnu2>c> alpha > beta >nn>xn>vn>tt>ct) 

function  - siaulate  error  occurrence  tine  and  correction  tine 

and  then  conpute  bawesian  estinate  of  nean 
correction  tine*  optinun  correction  linit  tiner 
and  its  nininun  cost  per  unit  tine?  and  also 
provide  optinun  sanple  size 

usage  - call  ndl2  (iseed>r»xnul»xnu2>c»alpha»beta>nn»xn>vn>tt>et) 

paraaeters  iseed  - (input.)  an  integer  value  in  the  exclusive  ranee 

( 1 >2147483647) . iseed  is  replaced  bv  a new  iseed 
to  be  used  in  subseouent  calls, 
r - (input.)  error  occurrence  rate 
xnul  - (input.)  mean  correction  tine  in  phase  1 

xnu2  - (input.)  nean  correction  tine  in  phase  2 

C - (input.)  a vector  of  length  3 

c(l)  - cost  Per  unit  tine  of  error  correction 

in  pase  1 

c(2)  - cost  Per  unit  tine  of  error  correction 
in  Pase  2 

c(3)  - sanpling  cost  per  sanple  size 
alpha  - (input.)  a vector  of  length  2 

alpha(l)  - shape  paraneter  of  inverted  ganna 
prior  for  nean  correction  tine 
in  pase  1 

alpha(2)  - shape  paraneter  of  inverted  ganna 
prior  for  error  occurrence  rate 
beta  - (input.)  a vector  of  length  2 

beta(l)  - scale  paraneter  of  inverted  ganna 
prior  for  nean  correction  tine 
in  pase  1 

beta (2)  - scale  paraneter  of  inverted  aanna 
prior  for  error  occurrence  rate 
nn  - (output.)  optinun  sanple  size 

xn  - (output.)  bavesian  estinate  of  ereor  occurence  rate 
vn  - (output.)  bavesian  estinate  of  nean  correction 
tine  in  pase  1 

tt  - (output.)  optinun  correction  linit  tine 
ct  - (output.)  nininun  cost  per  unit  tine  . 
read,  subroutines  - ggub>  optnn 


l 


subroutine  ndl2  ( iseed. r.xnul.xnu2.c.alpha>beta.nn.xn.wn.tt.ct> 


i 
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TABLE  C. 20  (Continued) 


dimension  c<3) >alpha(2) »beta(2) >xv(2) 

xx*©»0 

ww*>0»0 

do  5 nn=lflOO 

call  daub  (iseed>2rxy) 

xx=xx-r*aloa<xy<  1 ) ) 

vv*yy-x»u 1 * a 1 od < xy < 2 ) > 

if  <nn  .eo.  1)  ao  to  5 

xn=<xx+beta(2>  >/<alrha<2>+f loat<nn-l ) ) 

a=c<3)*f loat(nn) 

b*xx+xn+yy 

call  optnm  (nnryyrcralpha>betarxi»u2rarbrynrtlrcl) 

a«c(3>*float(nn+l) 

b=xx+2 . 0*xn+yy+yn 

call  optmih  <nn»yy »c»alpha»beta»x*u2» arbr ynrt2rc2) 
rrint  100»nn»xn»yn»tl »cl »t2»c2 
100  for»at<i5'6el5.5> 
if  <tl-t2)  11*22*5 
22  if  <cl-c2>  11*11*5 
5 continue 
11  tt=tl 
ct*cl 
return 
end 
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TABLE  C. 22  SUBROUTINE  GGUB 


subroutine  gsub  ( iseed.n* r )■ 


-fro*  imsl- 


function 

usage 

parameters  iseed 


basic  unifora  (0*1)  pseudo-random  number 
generator 

call  ggub  ( iseed *n*r) 

(input.)  an  integer  value  in  the  exclusive 
range  (1 *2147483647) . iseed  is  replaced  bv 
a new  iseed  to  be  used  in  subseauent  calls, 
(input.)  no.  of  deviates  to  be  generated 
output  vector  of  length  n»  containing  the 
deviates  in  (0*1) 


subroutine  ggub( iseed *n* r ) 
dimension  r(l) 

double  precision  z*dpm*dpn*dp 
data  dp**dpn/2147483<S47.dO*l.dO/ 
dp=dpn/ ( dpm+dpn ) 
z* iseed 
do  5 i-l*n 

z=dmod( 16807. dOCz.dpm* 

5 r(i)=z*dP 
iseed=z 
return 


TABLE  C.23 


SUBROUTINE  OPTMM 


e subroutine  ortat  (nn>vvicialrha»beta>xau2iaib>vn>itrct) 

C function  - compute  bauesian  estimate  of  acan  correction 

c timet  optimum  correction  limit  time*  and 

c its  minimum  cost  per  unit  time 

c usage  - call  optmm  (nn.uv.c* alpharbetarxmu2rarbrvnrttrct) 

c parameters  nn  - t input.)  no.  of  observations 

c w»  - (input.)  total  amount  of  observed  correction 

c time 

c c - ( input •>  a vector  of  length  3 

c C(l>  - cost  per  unit  time  of  error  correction 

c in  pase  1 

c c(2>  - cost  per  unit  time  of  error  correction 

c in  ease  2 

c c<3>  - sampling  cost  per  sample  size 

c alpha  - (input.)  a vector  of  length  2 

c alpha(l)  - shape  parameter  of  inverted  gamma 

c prior  for  mean  correction  time 

c in  pase  1 

c alpha<2)  - shape  parameter  of  inverted  gamma 

c prior  for  error  occurrence  rate 

c beta  - (input.)  a vector  of  length  2 

c beta(l)  - scale  parameter  of  inverted  gamma 

c prior  for  mean  correction  time 

c in  pase  1 

c beta(2)  - scale  parameter  of  inverted  gamma 

c prior  for  error  occurrence  rate 

c xmu2  - (input.)  mean  correction  time  in  pase  2 

c a - (input.)  constant  of  a in  eouation  (a-1) 

c b - (input.)  constant  of  b in  eouation  (a-1) 

c vn  - (output.)  bauesian  estimate  of  mean  correction 

c time  in  pase  1 

c tt  - (output.)  optimum  correction  limit  time 

c ct  - (output.)  minimum  cost  per  unit  time 


subroutine  optmm  (nn.uu. c» alpha. beta. xmu2tatb?wnt ttrrt) 

dimension  c<3> talpha(2) >beta(2) 

aa*c(2)*b-a 

bb-(c(l)*b-a)/xmu2 

gl«vv+beta(l ) 

aO*alPha( 1 )+f loat(nn-l ) 

un-gl/aO 

S2-bb/(aa+vn*c(2)-c(l>> 

g3-(a0+1.0)/gl 

t«(s3/s2-1.0)*sl 

if  (t  .It.  0.0)  go  to  22 

do  5 i*l » 20 

rt=g3/ ( 1 . 0+t/dl ) 

ggt-gl/a0*(l .0-(1.0+t/sl)**(-a0)> 


gt-(l .0+t/gl )**(-a0-l .0) 

xl=aa+(c(2)-c(l> >*agt 

f0»rt*xl+(c(2)-c( 1 > )*st-bb 

rrt«-g3/gl/( 1 .0+t/dl )**2 

ff»rrt*xl 

h—fO/ff 

t-t+h 

print  lOOtitt 
100  format(iS»elS.5) 


C-32 


TABLE  C. 23 


(Continued) 


If  <abs<h/t>  .1*.  0.00001)  so  to  11 
S continue 
11  tt«t 

et»(c<l>-c<2>*x»u2*rt>/< 1 . 0-x»u2*rt> 
return 

22  tt^o.O 

ct-<»+c<2)*x»u2>/<b+xeu2> 

return 


IBIS  PISE  IS  BBT  QUALITY  PRACTICABU 
rao«  oqpy  furbished  ro  doc  
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